clear; clc;
% 固定 A=0.35，f=0.72，扫描电场 E，绘制分岔图和 LE_E vs E

% 参数设定
a     = 0.7; 
b     = 0.8; 
c     = 0.1; 
alpha = 0.25; 
beta  = 0.01; 
delta = 0.1;

y0 = [0.2; 0; 0.01]; % 初始条件

% 固定外部激励参数：A和f
A_fixed = 0.35;
f_fixed = 0.72;

T_total = 500;   % 总模拟时间
T_trans = 100;    % 舍去瞬态时间
dt = 0.01;       % 用于 LE_E 估计的重新归一化间隔
h = 0.01;        % RK4 积分步长

% 扫描电场 E 的范围
E_range = linspace(0, 6, 300);

% 预分配存储数据
bifurcationE = [];   % 每个 E 下的 x_peak 数据，格式 [E, x_peak]
LE = zeros(size(E_range));

for i = 1:length(E_range)
    E_val = E_range(i);
    params = struct('a', a, 'b', b, 'c', c, 'alpha', alpha, ...
                    'beta', beta, 'delta', delta, 'A', A_fixed, ...
                    'f', f_fixed, 'E', E_val);
    
    % 使用修改后的 hybridNeuron_E 数值积分
    tspan = [0, T_total];
    [t, Y] = rk4(@(t,y) hybridNeuron_E(t, y, params), tspan, y0, h);
    
    % 舍去瞬态数据
    idx = t > T_trans;
    t_post = t(idx);
    x_post = Y(idx,1);
    
    % 根据外部激励周期采样：这里每个周期采样一次 x 的峰
    % 使用 findpeaks 检测 x(t) 的局部极大值
    [pks, ~] = findpeaks(x_post, t_post);
    
    % 将当前 E 值对应的所有峰值存入 bifurcationE
    bifurcationE = [bifurcationE; repmat(E_val, length(pks), 1), pks];
    
    % 计算当前 E 下的最大 Lyapunov 指数
    LE(i) = LE_E(params, y0, T_total, dt, h);
end

% (a) 分岔图：x_peak vs E
subplot(1,2,1);
plot(bifurcationE(:,1), bifurcationE(:,2), 'r.', 'MarkerSize', 5);
xlabel('E');
ylabel('x_{peak}');
title('Bifurcation diagram (f=0.72, A=0.35)');

% (b) LE vs E
subplot(1,2,2);
plot(E_range, LE, 'g-', 'LineWidth', 1.5); hold on;
yline(0, 'r--', 'LineWidth', 1.5);
xlabel('E');
ylabel('LE');
title('LE_E vs E (f=0.72, A=0.35)');
